Stationary spatial structures in reaction-subdiffusion 
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We discuss stationary concentrations of reactants in an A + B — > reaction under subdiffusion and 
show that they are described by stationary reaction-diffusion equations with a nonlinear diffusion 
term. We consider stationary profiles of reactants' concentrations and of reaction zones in a flat 
subdiffusive medium fed by reactants of both types on its both sides (a subdiffusive gel reactor) . 
The behavior of the concentrations and of the reaction intensity in subdiffusion differs strikingly 
from those in simple diffusion. The most important differences correspond to the existence of 
accumulation and depletion zones close to the boundaries and to non-monotonous behavior of the 
reaction intensity with respect to the strength of the minor source. The implications of these results 
for other situations are also discussed. 

PACS numbers: 05.40.Fb, 82.33.Ln 



Many phenomena in systems out of equilibrium are 
described as reactions between diffusing species. Apart 
from chemistry, examples are trapping and annihila- 
tion of excitons or recombination of charge carriers in 
physics, or predator-prey relations in ecology. For nor- 
mal diffusion the situation is described in terms of 
reaction-diffusion equations. Many systems however ex- 
hibit anomalous diffusion with a mean square displace- 
ment scaling as (r 2 (t)') oc t a , with < a < 1 in the 
subdiffusive case and a > 1 in the superdiffusive case, 
which is not described by a diffusion equation [J Q. 
Subdiffusion can often be modeled within the frame- 
work of continuous time random walks (CTRW) with a 
heavy-tailed waiting time density function ip(t) oc t~ x ~ a , 
which yields a fractional subdiffusion equation. In anal- 
ogy to reaction-diffusion, fractional reaction-subdiffusion 
equations have been proposed, where an additional frac- 
tional time derivative acts either on the spacial Lapla- 
cian [H, Q| or on both the spatial Laplacian and the re- 
action term [H, 0]. Refs.0, 0] showed that the reaction- 
subdiffusion equations do not follow by simply changing 
a diffusion operator in a reaction-diffusion equation for 
a subdiffusion one. The situation here was pertinent to 
initial-condition problems. A rather general approach to 
reaction-subdiffusion equations was proposed in Ref.0], 
where, after deriving the set of reaction-subdiffusion 
equations the authors analyzed the properties of Tur- 
ing instability in corresponding reaction-subdiffusion sys- 
tems. In what follows we discuss the simplest boundary 
condition problem in which nontrivial stationary spacial 
structures appear, namely the irreversible A + B — > 
reaction [f|, and focus on stationary concentration pro- 
files under given, steady concentrations of reactants on 
the boundaries (Ref.[5[ started with one reactant, A or 
B, on each side). We show that the spacial structures 
under reaction-subdiffusion are strikingly different from 
ones emerging in reaction-diffusion. We first give the 
derivation of subdiffusion-reaction equations and station- 
ary subdiffusion-reaction equations by generalizing the 



scheme put forward in Refis.0,|||. As in 0, subdiffusion 
is considered within a CTRW scheme, and the reaction 
locally follows the mass action law. We then proceed by 
discussion of stationary forms of concentrations' distri- 
butions and of those of reaction zones. 

In a CTRW a particle arriving at a site i at time t' stays 
there for a sojourn time t, which is given by the probabil- 
ity density function ip(t). Leaving the site it makes a step 
with probability 1/2 in either direction. In the following 
discussion we confine ourselves to the one-dimensional 
situation, however the generalization to whatever other 
geometry is quite evident. 

The generalized reaction-diffusion equations are based 
on two balance conditions. The balance equation for A- 
particles at each site reads: 

A i (t)=j?(t)-jr(t) + R i {A,B} (1) 
= \h-itt) + 5^1 (*) - hit) + RdA, B}. (2) 

where j~(t) is the loss flux of A-particles at site i, i.e. 
the probability for an A-particle to leave i per unit time, 
jf{t) is the gain flux at the site, and Ri{A, B} — —kAiBi 
is a reaction term, describing particles' loss due to reac- 
tion. Since in our case the equations for A and B particles 
are symmetric, we concentrate on the equations for the 
A particles. A generalized reaction-diffusion equation is 
a combination of the continuity equation, Eq.([2|), and 
the equation for the loss fluxes j~(t) following from the 
assumption about the distribution of sojourn times tp(t) 
and survival probability P(i, to). 

According to the sojourn time distribution, the loss 
current for site i at time t is connected to the gain cur- 
rent for the site at all previous times and with the sur- 
vival probability. Namely, the particles which leave site 
i at time t (making a step from i to one of its neigh- 
bors) were either at site i from the very beginning (and 
survived), or arrived there at some time t' < t and sur- 
vived until t. The probability density to make a step 
at time t, having arrived at t' , is given by the wait- 
ing time distribution ip(t — t'). We have then j~(t) = 



2 



il>(t)Pi(t, 0)A i (0) +J^(t-t')P i {t,t')jf (t')dt', which, by 
using Eq.|T|), can be rewritten in the form 

j-(t) = i/>(]t)Pi(t,o)Mo) + I ^(*- O«(*.*0 x 

Jo 

x AiW+jiW + kAi^Biit'^dt'. (3) 



The survival probability of A at i is given by the classical 
kinetic equation j%Pi{t) = —kBi(t)Pi(t) and depends on 
the time-dependent B-concentration via 



Pi(t,ta) = exp -k / Bi(t')dt' 



(4) 



At this stage we also can assume the concentrations to 
be slowly changing is space, and change to a continuous 
coordinate x = ai, with a being the lattice spacing: 

a 2 

A[x,t) = —Aj-(x,t) - kA(x,t)B(x,t). (5) 

Eq.([5]) together with Eqs.([3]) and ((4]) and their coun- 
terparts for B-concentration give us the full system of 
equations for time-dependent concentrations. This sys- 
tem can be transformed to a special case of reaction- 
subdiffusion equations considered in Ref. 0. However, 
at this stage it is easier to proceed with the equations in 
our form. 

In our previous discussion we considered an initial con- 
dition problem, where A-particles were introduced at 
time t = into the system and followed the evolution 
of their concentration. In the case of normal diffusion 
the situation in steady state (achieved, say, if the con- 
centration of the particles at boundaries of the system is 
fixed by external sources) is given by the same reaction- 
diffusion equations, with the time derivatives at the left 
hand side put to zero. For the case of subdiffusion the 
situation is a bit more involved. 

Let us assume that in the course of time the system 
achieves a steady state characterized by constant con- 
centrations A(x) and B{x). Such a steady state is main- 
tained through the particles' sources at the boundaries of 
the system, no particle sources exist in the interior. Let 
us label the particles according to the time to they were 
introduced into the system, so that e.g. A(x,t\to)dto is 
the concentration at point x at time t of A-particles in- 
troduced between to and to + dto (a partial concentration 
of A). The partial concentration A(x, to|io) of newly in- 
troduced particles is zero everywhere in the interior of 
the system. The overall concentration of A-particlcs at 
site x is given by the integral 



A(x) 



A(x,t\t )dt . 



(6) 



In a steady state A(x,t\to) can only be a function of the 
difference of the time-arguments, i.e. of the elapsed time 



t e = t — to so that A(x, t\to) — A(x, t — to) and the equa- 
tions for the partial concentrations A(x,t\to) are given 
by Eqs.([5|) with t changed to t e . Moreover, the overall 
concentration A(x) is given by A(x) = J Q °° A(x,t e )dt e . 
Since A(x) and B(x) are time- independent, the sur- 
vival probabilities Pa(x, t, t') and Pb(x, t, t') in reaction- 
subdiffusion equations are the functions of the differ- 
ences of their time-arguments so that PA{x,t,t') = 
exp [— kB(x)(t — t')). The integral in the equation for 
the flux now takes the form of a convolution: 



j (x, t\t Q ) = ip(t)P A (x, t - t )A(x, t \t ) 
+ / ip(t-t')P A (x,t-t') X 



(7) 



x A(x,t'\to)+j~{x,t'\to) + kA(x,t'\t )B(x) dt' 

where j~(x, t'\to) are the loss fluxes for those A-particles 
which were introduced into the system at time to- 
We now pass to the Laplace domain with respect 
to t e and denote A(x,u) = J t A{x, t\to) exp[— u[t — 
to)]dt. The Laplace transform of the product ^>(t,x) = 
ip(t)exp[—kB(x)t] is given by the shift theorem and is 
equal to ^(u, x) — tp(u + kB(x)), so that 

~ = [u + kB(xM(u + kB(x)) A 

' 1 - ^{u + kB(x)) 

Inserting this into the first equation we get in the Laplace 
domain 



uA(x, u) — A(x, to\to) = 
a 2 [u + kB(x)]^{u + kB(x)) ~ 



(9) 



2 1 - ij)(u + kB{x)) 



A{x 1 u) — kA(x, u)B(x), 



where A(x, to \tq) differs from zero only at the boundaries. 
With A(x) = A(x,0) the stationary concentration A{x) 
in the interior of the system is given by 



a 2 kB{x)i){kB{x)) 



1 - ${kB(x)) 



A{x) - kA{x)B(x) = 0, (10) 



with the boundary conditions corresponding to the 
given concentrations on the boundaries. For a Marko- 
vian case of regular diffusion, corresponding to ip(t) = 
r _1 exp(— t/r), one has i/j(u) = 1/(1 + ut) so that this 
equation reduces to (a 2 /2t)AA(x) — kA(x)B(x) = 0, a 
usual stationary reaction-diffusion equation. In the non- 
Markovian case, corresponding to subdiffusion, the wait- 
ing time distribution in the Laplace domain is ip(u) ~ 
1 — (T?i)°T(l — a) for small u and Eq. (fT0|) reads 

° 2 1 

AB{xf- a A{x) ~ k a A(x)B(x) = 0. (11) 



2 r Q r(l - a) 



It is easy to see that the Markovian equation is a spe- 
cial case of Eq. (TTTj) for a = 1. The combination D a — 



3 



a 2 /2r Q r(l — a) stands for a (generalized) diffusion coef- 
ficient. The full system of steady state equations is given 
by Eq.(fTU|) and the corresponding equation for B. It is 
interesting to stress that the system of equations with 
additional temporal operator acting on the Laplacian in 
the case of initial-condition problem turns to a system of 
nonlinear reaction-diffusion equations with nonlinearity 
in the diffusion term for a stationary boundary-condition 
problem. 

As an example we consider a system on an interval 
(0, 1) with given reactants' concentrations on the bound- 
aries of the interval. The physical system we have in 
mind is a gel reactor or a porous medium in contact with 
two well-mixed reservoirs on both sides. The concentra- 
tions A(Q) — B(l) — 1 at the major sources are fixed. 
For the sake of simplicity we consider here symmetric 
situation with A(Q) = B(l) and A(l) = B(0). Due 



to symmetry the ^-concentration is then always sym- 
metric to the ^-concentration, i.e. B(x) = A(l — x). 
The concentration B(0) = A(l) will be called the mi- 
nor source strengths (except for the symmetric case 
A(l) = B(0) = 1). In Fig.l we show numerical results 
for the steady-state Eqs. tfTTj) obtained by semi-implicit 
relaxation algorithm. The results for subdiffusion (here 
with moderate a = 0.9) are compared to the ones for 
normal diffusion (a = 1 and hence linear diffusion opera- 
tor). Shown is the behavior of the A-concentrations and 
the one of the reaction intensity kA(x)B(x). The param- 
eters are: k = 0.01, D a = l/2r(l - a) for a = 0.9 and 
D a = 1/2 for a = 1. The A(0) = B(l) concentration is 
fixed to be A(0) = 1, the other concentration varies from 
5(0) = A{1) = 1 (symmetric case, when the CTRW- 
reactor separates two stoichiometric reacting mixtures) 
to B(0) = A(l) = 10~ 4 . 
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FIG. 1: Stationary particle concentration profile A(x) for 
a = 1 (upper panel) and for a = 0.9 (lower panel); A(l) = 1 
(dashed), 0.5 (dash-dot-dot), 5-10" 2 (dash-dot), 5-10" 3 (dot- 
ted), and 1 • 10 -4 (solid line), see text for details. Note the 
growth of the peak and its subsequent decay and outwards 
motion with decreasing minor source strength for the subdif- 
fusive case; the behavior of the particle concentration profile 
in regular diffusion is monotonous with respect to the strength 
of the minor source and tends to a limiting form for A(l) — > 0. 



FIG. 2: Stationary reaction intensity profile kA(x)B(x) for 
a = 1 (upper panel) and for a — 0.9 (lower panel). The pa- 
rameters are the same as in Fig.l, namely: A(l) = 1 (dashed 
line) , 0.5 (dash-dot-dot), 5 ■ 10" 2 (dash-dot), 5 • 10" 3 (dot- 
ted), and 1 ■ 10~ 4 (solid line). In the diffusive case the sta- 
tionary reaction intensity profile approaches a limiting form 
under decreasing minor source strength. For subdiffusion- 
reaction the height of the profile shows nonmonotonous be- 
havior with respect to this parameter. Note the difference in 
vertical scales! 
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In the symmetric case both reaction-diffusion and 
reaction-subdiffusion situations correspond to a very sim- 
ilar behavior with maximal concentrations in the regions 
close to the boundaries where the system is fed by reac- 
tants. For asymmetric boundary conditions however the 
behaviors of the concentrations in regular diffusion and 
in subdiffusion differ strongly. One of the most marked 
differences corresponds to a strong non-monotonicity of 
the concentration close to the major source (side with the 
higher concentration of A) indicating for accumulation of 
A particles in the interior of the subdiffusive medium. Its 
counterpart on the other side of the system is a depletion 
zone (corresponding to the symmetric accumulation zone 
for B). The peak and the depletion zone are much more 
pronounced for smaller a; the choice of moderate a = 0.9 
was caused by our whish to show the behavior pertinent 
to diffusion and to subdiffusion on the same scale. It 
is important to note that the dependence of the height 
of the accumulation peak on the strength of the minor 
source is nonmonotonous: The reduction of the minor 
source strength A(l) = B(0) leads first to its growth, 
and then to its outwards motion accompanied by decay. 

Let us now turn to reaction intensities, Fig. 2. For 
the symmetric case the reaction takes place mostly close 
to the boundaries of the system. For smaller A(l) the 
reaction zone starts to form in the middle of the sys- 
tem. However, also here striking differences between the 
reaction-diffusion and the reaction-subdiffusion cases are 
seen. In the reaction-diffusion case the dependence of the 
form of the reaction zone on A(l) is weak for small A(l), 
and there exists a clear limiting form for A(l) = 0. This 
behavior is known and is used in the time-scale separa- 
tion approach of Refs.fiol. 11 1 based on the quasistatic 



approximation. For reaction-subdiffusion the behavior 
of the reaction zone with respect to its height is non- 
monotonous. When lowering A(l), the maximum of re- 
action intensity first gets higher and then starts to lower; 
and the distribution as a whole broadens. The reason for 
this is quite evident. The stationary reaction zone exists 
only if it is fed by A- and B-reactants on the correspond- 
ing sides. Both in the diffusion and in the subdiffusion 
case the reaction zone is the higher and the narrower the 
larger is the particles' inflow into the reaction area. This 
inflow is governed by the effective diffusion coefficient 
of the corresponding reactants, which, in the subdiffu- 
sive case, depends on the concentration of the reacting 
counterpart. Since the effective mobility of subdiffusing 
species decays in the course of time (the number of steps 
per unit time goes as the effective diffusion in sub- 

diffusion is caused by the reaction itself: it corresponds to 
a reaction-induced diffusion term, just as in the reactions 
between immobile species , [HI]. For zero A(l) = B(0) 
concentration the effective diffusion coefficient vanishes 



on the corresponding side of the system preventing the 
inflow of reactants from their major sources into the in- 
terior of the system. The reaction zone blurs and fades 
out. In this case no stationary front exists. This effect is 
also clearly seen when considering the time evolution of 
concentrations which can be done by discussing the prop- 
erties of the inverse Laplace transform of Eq.lfTOll. This 
means that the adiabatic approximation of Refs. 

[13, EI 

fails in subdiffusion, and the analysis of the front's mo- 
tion in this case has to be done anew. 

Let us summarize our findings. We discussed the sta- 
tionary form of reactants' concentrations and of reac- 
tion zones in the A + B — > reaction in a subdiffu- 
sive medium fed by reactants on its both sides. We 
show that the behavior of the concentration and of the 
reaction intensity profiles in subdiffusion differs strik- 
ingly from those in simple diffusion. The most important 
differences correspond to the existence of accumulation 
and depletion zones close to the boundaries and to non- 
monotonous behavior of the reaction intensity with re- 
spect to the strength of the minor source. These results 
have implications for the time-dependent front forms in 
A + B — > reaction, since the adiabatic approxima- 
tion used in the reaction-diffusion case would fail in the 
reaction-subdiffusion one. 
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